Add Folders to Path

%%time
import sys, os

# get current directory
path = os.getcwd()

# get parent directory
parent_directory = os.path.sep.join(path.split(os.path.sep)[:-4])

# add utils folder to current working path
sys.path.append(parent_directory+"/subfunctions/utils")

# add integration folder to current working path
sys.path.append(parent_directory+"/subfunctions/integration")

# add HyperbolicsOECS folder to current working path
sys.path.append(parent_directory+"/demos/AdvectiveBarriers/HyperbolicOECS")

Wall time: 0 ns

Overview

In the following notebok we extract elliptic OECS in the Agulhas region from the AVISO dataset. The notebook is structured as follows:

1. Import data from the file 'AVISO.mat' stored in the folder 'data'.
   

2. Define computational parameters (such as the number of cores) and data.
   

3. Define spatio-temporal domain.
   

4. Interpolate velocity from the (discrete) gridded data.
   

5. Hyperbolic OECS:

   • Compute rate of strain tensor $ \mathbf{S}(\mathbf{x}, t) $ over meshgrid.

   • Compute eigenvalues $ s_i(\mathbf{x},t) $ and eigenvectors $ e_i(\mathbf{x},t) $ (i = {1,2}) of rate of strain tensor.

   • Compute the sets of isolated local maxima of $ \mathbf{S}(\mathbf{x}, t) $ (=objective saddle-points)

   • Compute repelling OECS as tensorlines tangent to the eigenvectors $ \mathbf{e_1}(\mathbf{x}) $ of the rate of strain tensor. The tensorlines are launched from isolated local maxima (objective saddle-points). \begin{equation} \dfrac{d}{ds}\mathbf{e}(\mathbf{x}) = \mathbf{e_1}(\mathbf{x}) \end{equation}

   • Compute attracting OECS as tensorlines tangent to the eigenvectors $ \mathbf{e_2}(\mathbf{x}) $ of the rate of strain \begin{equation} \dfrac{d}{ds}\mathbf{e}(\mathbf{x}) = \mathbf{e_2}(\mathbf{x}) \end{equation}

   • Stop the integration of the tensorlines as soon as $ |s_2(\mathbf{x},t)| $ stops to be monotonically decreasing or $ |s_2(\mathbf{x},t)| < d_{hyperbolicity} $, where $ d_{hyperbolicity} $ is a threshold imposed on the rate of attraction/repulsion (=hyperbolicity) of the attracting/repelling (=hyperbolic) OECS.

Import Data

%%time
import scipy.io as sio

#Import velocity data from file in data-folder
mat_file = sio.loadmat('../../../../data/Aviso/AVISO.mat')

U = mat_file['u']
V = mat_file['v']
x = mat_file['x']
y = mat_file['y']
time_data = mat_file['t']

Wall time: 116 ms

Computational parameters and data

Here we define the computational parameters and the data.

import numpy as np

# number of cores to be used for parallel computing
Ncores = 8

# time resolution of data
dt_data = time_data[0, 1]-time_data[0,0]

# periodic boundary conditions
periodic_x = False
periodic_y = False
periodic = [periodic_x, periodic_y]

# unsteady velocity field
bool_unsteady = True

# defined domain
defined_domain = np.isfinite(U[:,:,0]).astype(int)

# set nan values

## compute meshgrid of dataset
X, Y = np.meshgrid(x, y)

## resolution of meshgrid
dx_data = X[0,1]-X[0,0]
dy_data = Y[1,0]-Y[0,0]

delta = [dx_data, dy_data]

Spatio-temporal domain

Here we define the spatio-temporal domain over which to consider the dynamical system.

%%time
# Time
t_OECS = 0

# store time in array
time = np.array([t_OECS])

# longitudinal and latitudinal boundaries (in degrees)
xmin = -3
xmax = 1
ymin = -32
ymax = -24

# make sure that domain is in the data
assert (xmax <= np.max(X) and xmin >= np.min(X) and ymin >= np.min(Y) and ymax <= np.max(Y) and np.min(time_data) <= t_OECS <= np.max(time_data)),"The domains you are chooising are outside the domain of the data!!!!! --> redefine spatial/temporal domain"
    
# spacing of meshgrid (in degrees)
dx = 0.025
dy = 0.025

x_domain = np.arange(xmin, xmax + dx, dx)
y_domain = np.arange(ymin, ymax + dy, dy)

X_domain, Y_domain = np.meshgrid(x_domain, y_domain)

Wall time: 1.03 ms

Interpolate Velocity

In order to evaluate the velocity field at arbitrary locations and times, we must interpolate the discrete velocity data. The interpolation with respect to time is always linear. The interpolation with respect to space can be chosen to be "cubic" or "linear". In order to favour a smooth velocity field, we interpolate the velocity field in space using a cubic interpolant.

%%time
# Import interpolation function for unsteady flow field
from ipynb.fs.defs.Interpolant import interpolant_unsteady

# set nan values to zero so that we can apply interpolant. Interpolant does not work if the array contains nan values
U[np.isnan(U)] = 0
V[np.isnan(V)] = 0

# Interpolate velocity data using cubic spatial interpolation
Interpolant = interpolant_unsteady(X, Y, U, V, time_data, method = "cubic")

Wall time: 87.8 ms

Rate of strain tensor over meshgrid of initial conditions

The rate of strain tensor $ S(\mathbf{x}, t) $ at time $ t $ is computed by iterating over meshgrid. The rate of strain tensor at point $ \mathbf{x} $ at time $ t $ is computed from the gradient of the velocity field by using an auxiliary meshgrid. 'aux_grid' specifies the ratio between the auxiliary grid and the original meshgrid. This parameter is generally chosen to be between $ [\dfrac{1}{5}, \dfrac{1}{10}] $.

%%time
# Import package for progress bar
from tqdm.notebook import tqdm

# Import package for parallel computing
from joblib import Parallel, delayed

# Import gradient of velocity function
from ipynb.fs.defs.gradient_velocity import gradient_velocity

# Import package which checks particle location
from ipynb.fs.defs.check_location import check_location

# Import Rate of Strain function
from ipynb.fs.defs.RateStrain import RateStrain

# Import eigenvalue/eigenvector calculator
from ipynb.fs.defs.eigen import eigen

# Define ratio of auxiliary grid spacing vs original grid_spacing
aux_grid_ratio = .2 # [1/5, 1/10]
aux_grid = [aux_grid_ratio*(X_domain[0, 1]-X_domain[0, 0]), aux_grid_ratio*(Y_domain[1, 0]-Y_domain[0, 0])]

def parallel_S(i):
    
    S_parallel = np.zeros((X_domain.shape[1], 2, 2))*np.nan
    
    for j in range(S_parallel.shape[0]):
        
        x = np.array([X_domain[i,j], Y_domain[i,j]])
        
        # only compute rate of Strain for trajectories starting region where velocity field is defined
        if check_location(X, Y, defined_domain, x)[0] == "IN":
        
            # Compute gradient of velocity
            grad_vel = gradient_velocity(t_OECS, x, X, Y, Interpolant, periodic, defined_domain, bool_unsteady, dt_data, delta, aux_grid)
        
            # Compute rate of strain at 'x'
            S_parallel[j, :, :] = RateStrain(grad_vel)
    
    return S_parallel

S = np.array(Parallel(n_jobs=Ncores, verbose = 0)(delayed(parallel_S)(i) for i in tqdm(range(X_domain.shape[0]))))

Wall time: 11.3 s

Eigenvalues/Eigenvectors of rate of strain tensor

We now compute the properties of the rate of strain tensor 'S' such as the eigenvalues 's1', 's2' and eigenvectors 'eigenv1', 'eigenv2'.

# Import eigenvalues/eigenvectors calculator
from ipynb.fs.defs.eigen import eigen

eig1 = S[:,:,0,0].copy()*np.nan
eig2 = S[:,:,0,0].copy()*np.nan
e1 = np.zeros((S.shape[0], S.shape[1], 2))
e2 = np.zeros((S.shape[0], S.shape[1], 2))

#iterate over meshgrid
for i in range(X_domain.shape[0]):
    
    for j in range(Y_domain.shape[1]):
            
        x = [X_domain[i,j], Y_domain[i, j]]
        
        # compute eigenvalues of rate of strain
        eig1[i,j], eig2[i,j], e1[i,j,:], e2[i,j,:] = eigen(S[i,j,:,:])

Find objective Saddle-points

Objective saddles coincide with local maxima in the maximum eigenvalue field $ s $ of the rate of strain tensor.

# import function to compute local maxima
from ipynb.fs.defs.loc_max import _loc_max

# minimum distance between objective saddles (=local maxima in the s2-field)
min_distance = .75

# indices, positions and value (of s2-field) of the objective saddles
loc_idx_x, loc_idx_y, loc_max_x, loc_max_y, loc_max_field = _loc_max(min_distance, X_domain, Y_domain, eig2)

Shrinklines (Repelling OECS)

Repelling LCS can be sought among trajectories of the differential equation:

\begin{equation} \mathbf{x}'_0(s) = \mathbf{e}_1(\mathbf{x}_0;t), \label{eq: shrinklines} \end{equation}

with $ \mathbf{e}_1 $ denoting the eigenvector associated to the weakest eigenvalue $ s_1 $ of $ S(\mathbf{x},t) $. The non-orientable vector field is well defiend away from tensorline singularites (points where $ s_1 = s_2 $). The most repelling shrinklines mark initial positions of repelling OECSs. Repelling OECSs can therefore be located as trajectories of eq. \ref{eq: shrinklines} that have locally the largest averaged $ s_2(\mathbf{x},t) $ among all neighbouring shrinklines.

from ipynb.fs.defs.tensorlines import _tensorlines

# step-size used for integration with respect to parameter 's'
step_size = 0.01

# threshold distance to locate local maxima in the 'eig2'-field
max_distance = 0.75

# maximum length of shrinkline
max_length = 3

# Minimum threshold on rate of repulsion of shrinkline 
hyperbolicity = .15

shrinklines = _tensorlines(X_domain, Y_domain, eig2, e1, max_distance, max_length, step_size, hyperbolicity)

Stretchlines (Attracting OECS)

Attracting LCS can be sought among trajectories of the differential equation:

\begin{equation} \mathbf{x}'_0 = \mathbf{e}_2(\mathbf{x}_0;t), \label{eq: stretchlines} \end{equation}

with $ \mathbf{e}_2 $ denoting the eigenvector associated to the strongest eigenvalue $ s_2 $ of $ S(\mathbf{x},t) $. The non-orientable vector field is well defiend away from tensorline singularites (points where $ s_1 = s_2 $). The most attracting stretchlines mark initial positions of attracting OECSs. Attracting OECSs can therefore be located as trajectories of eq. \ref{eq: stretchlines} that have locally the smallest averaged $ s_1(\mathbf{x},t) $ among all neighbouring stretchlines.

from ipynb.fs.defs.tensorlines import _tensorlines

# step-size used for integration with respect to parameter 's'
step_size = 0.01

# threshold distance to locate local minima in the 's_1' field of 'S'
max_distance = 0.75

# maximum length of stretchline
max_length = 5

# Minimum threshold on rate of attraction of stretchline 
hyperbolicity = 0.15

stretchlines = _tensorlines(X_domain, Y_domain, eig2, e2, max_distance, max_length, step_size, hyperbolicity)

################################ PLOT HYPERBOLIC OECS ################################

import matplotlib.pyplot as plt

# generate figure and axis object
fig = plt.figure(figsize = (16, 12), dpi = 600)
ax = plt.axes()

# Plot eig2-field
cax = ax.contourf(X_domain, Y_domain, eig2, levels = 600, cmap = "gist_rainbow_r")

cbar = plt.colorbar(cax)
cbar.set_ticks(np.arange(0, 1.5, 0.1))
cbar.set_label(r'$ |s_1| $', rotation = 0, labelpad = 20, fontsize = 16)

# Plot attracting OECS  
for i in range(len(stretchlines[0])):
    ax_attracting = ax.plot(stretchlines[0][i], stretchlines[1][i], c = 'k', linewidth = 2.5)

# Plot repelling OECS
for j in range(len(shrinklines[0])):
    ax_repelling = ax.plot(shrinklines[0][j], shrinklines[1][j], c = "w", linewidth = 2.5)

# set limits
ax.set_xlim(np.min(X_domain), np.max(X_domain))
ax.set_ylim(np.min(Y_domain), np.max(Y_domain))

# Set title
ax.set_title('Repelling/Attracting OECSs', fontsize = 20)

# Plot legend
import matplotlib.patches as mpatches
attracting = mpatches.Patch(color='k', label='attracting OECS')
repelling = mpatches.Patch(color='w', label='repelling OECS')
plt.legend(handles=[attracting, repelling], fontsize = 16,  loc = "lower left")
 
# Show plot
plt.show()

Advect attracting OECS

Attracting/repelling OECS represent short term attractor/repellors of particles in the flow field. They act as the eulerian counterpart to hyperbolic LCS. Hyperbolic OECS act as a short-term unstable/stable manifold of a saddle point. In order to highlight the validity of the extracted hyperbolic OECS and highlight the predictive power of attracting OECS for short time-intervals, we advect attracting OECSs over 7 days and study their effect on nearby particles. We seed around every objective saddle a circular blob of particle which also gets advected by the flow. As shown in the below figure, the attracting OECS play a crucial role in shaping the deformation of the blobs as they stretch and attract the blob originally centered around the objective saddle.

A comparison with the advected attracting OECS computed with the FastTensorlineComputation reveals identical structures.

# Import package for computing trajectories
from ipynb.fs.defs.integration_dFdt import integration_dFdt

# efine time-horizon over which to advect
time_advect = np.linspace(t_OECS, t_OECS+7, 50, endpoint = True)

# Define figure/axes
fig = plt.figure(figsize = (16, 10), dpi = 600)
ax = plt.axes()

# Iterate over all attracting OECS and plot
for i in tqdm(range(len(stretchlines[0]))):
    
    x_advected, y_advected = [], []
                    
    for j in range(len(stretchlines[0][i])):
        
        x = np.array([stretchlines[0][i][j], stretchlines[1][i][j]])
                    
        # Advect OECS over 7 days
        Fmap = integration_dFdt(time_advect, x, X, Y, Interpolant, periodic, defined_domain, bool_unsteady, dt_data, delta)[0]

        # extract end-point of advection
        x_end = Fmap[:,-1]
                    
        x_advected.append(x_end[0])
        y_advected.append(x_end[1])
    
    ax.plot(x_advected, y_advected, linewidth = 5, c = "r")
    
    # create circle around objective saddle and seed it with particle
    r = 0.25
    x_circle, y_circle = [], []
    for phi in np.linspace(0, 2*np.pi, 100):
        x_circle.append(loc_max_x[i]+r*np.cos(phi))
        y_circle.append(loc_max_y[i]+r*np.sin(phi))
    
    x_circle_advected, y_circle_advected = [], []
    # advect particles around objective saddle over 7 days
    for t in range(len(x_circle)):
        
        x = np.array([x_circle[t], y_circle[t]])
                    
        # Advect OECS over 7 days
        Fmap = integration_dFdt(time_advect, x, X, Y, Interpolant, periodic, defined_domain, bool_unsteady, dt_data, delta)[0]

        # extract end-point of advection
        x_end = Fmap[:,-1]
                    
        x_circle_advected.append(x_end[0])
        y_circle_advected.append(x_end[1])
        
    ax.fill(x_circle_advected, y_circle_advected, c = "b", alpha = 1)
        
# Set axis limits
#ax.set_xlim([-5, 1])
#ax.set_ylim([-32, -24])
                
# Set axis labels
ax.set_xlabel("long (deg)", fontsize = 16)
ax.set_ylabel("lat (deg)", fontsize = 16)

# Title
ax.set_title("Advected attracting OECSs", fontsize = 20)
plt.show();

References

[1] Serra, M., & Haller, G. (2016). Objective Eulerian coherent structures. Chaos: An Interdisciplinary Journal of Nonlinear Science, 26(5), 053110.